Mesoscopic supersolid of dipoles in a trap 
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A mesoscopic system of indirect dipolar bosons trapped by a harmonic potential is considered. 
The system has a number of physical realizations including dipole excitons, atoms with large dipolar 
moment, polar molecules, Rydberg atoms in inhomogenious electric field. We carry out a diffusion 
Monte Carlo simulation to define the quantum properties of a two-dimensional system of trapped 
dipoles at zero temperature. In dimensionless units the system is described by two control param- 
eters, namely the number of particles and the strength of the interparticle interaction. We have 
shown that when the interparticle interaction is strong enough a mesoscopic crystal is formed. As 
the strength of interactions is decreased a multi-stage melting takes place. Off-diagonal order in 
the system is tested using natural orbitals analysis. We have found that the system might be Bose- 
condensed even in the case of strong interparticle interactions. There is a set of parameters for 
which a spatially ordered structure is formed while simultaneously the fraction of Bose condensed 
particles is non zero. This might be considered as a realization of a mesoscopic supersolid. 



I. INTRODUCTION 

Bose-Einstein condensation (BEC) is a phenomenon 
that has been attracting great attention for a long time 
since its prediction in 1924 [l|, [2j]. Theoretical description 
of condensate properties is commonly based on mean- 
field Gross-Pitaevskii equation Q which applies to a di- 
lute Bose gas. This assumption is not fulfilled in systems 
with strong correlations. 

This is the case of a two-dimensional rather dense 
dipole system of indirect excitons in coupled quantum 
wells or in single quantum well in strong electric field [ij . 
At the densities much smaller than a~ 2 (with a being the 
characteristic excitons size) exchange effects are greatly 
suppressed by the dipole-dipole repulsion, so dipole exci- 
tons can be treated as Bose particles. The exciton con- 
finement can be created by inhomogeneous electric field 
or inhomogeneous deformation of semiconductor. The 
inhomogeneous electric field can be generated, for exam- 
ple, by a tip of scanning probe-microscope or by a profiled 
controlling gate (see [5] and references therein). 

The dipolar interactions are also important in atomic 
gases with large dipolar moment. The Bose-Einstein 
condensation has been achieved in 52 Cr atoms possessing 
large permanent dipolar moment [6|. In this system 
there is a competition between short-range interactions 
and long-range dipolar forces. Strong dipolar effects 
has been experimentally observed in trapped chromium 
atoms, see Ref. Large dipolar moment can be 

induced in polar molecules by applying an external 
electric field. Recently the quantum regime of 40 K 87 Rb 
fermionic polar molecules has been reached and dipolar 
collisions in such systems has been experimentally 
studied [III [H| . Rydberg atom (atom with one electron 
excited to a very high principal quantum number (l3j ) 
has a very large size and a very large polarizability, 
i.e. large dipole moment can be inducted in moderate 
electric field. An electric field align dipolar moments so 



that the interaction in 2D system has dipolar 1 /r 3 form. 
An optical trapping can be used to confine the system 
of Rydberg atoms [14| . 

A close relation between two quantum phenomena, 
BEC and superfluidity, suggests that a system might be 
supcrfluid if a part of it is Bose condensed. Moreover, 
at the same time the system can possess a crystal order, 
i.e. can be a supersolid. A finite superfluid signal was 
reported in a number of recent experiments with cold 
helium in solids (see [tsU^ ) . It is probable that the su- 
perfluid signal observed in experiments is due to presence 
of defects. We note that mesoscopic trapped systems are 
good candidates for being supersolid as no perfect com- 
mensurate crystalline structure can be formed and de- 
fects are intrinsically present in the system. We perform 
numerical simulations of a trapped system of excitons 
with dipolar interaction, to obtain structure properties 
and to study how the crystal-like order is formed in the 
system. 

In the regime when dipole-dipole repulsions are weak, 
analytical approaches based on mean field approxima- 
tion, provide a good description of the system, but fail 
when the density is large. Instead Monte Carlo methods 
have no such limitations and can be successfully applied 
to strongly interacting systems. In Ref. [23| properties 
of a classical dipolar system in a harmonic trap at low 
finite temperature were studied. Path Integral Monte 
Carlo (PIMC) simulation of trapped clusters were per- 
formed in works of Lozovik et al. [24| and Pupillo et 
al. [IH for dipolar interaction and for Coulomb inter- 
action by Filinov et al. [25|. PIMC technique is useful 
for simulation of the properties of quantum systems at 
a finite temperature, while the smaller is the tempera- 
ture the more difficult it is to obtain accurate results, 
which makes the simulation of the ground state of the 
system extremely difficult. Instead diffusion Monte Carlo 
method is well suited for studying ground state proper- 
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ties. In references (2ril - |29| two-dimensional (2D) systems 
of dipoles in the absence of an external potential were 
studied. It was shown that at large density a crystal 
is formed. Non-zero superfluid and condensate fractions 
have been found in finite-size crystals and in crystals with 
vacancies. 

In the present work we study effects of a harmonic 
confinement in a 2D geometry. Our motivation for ex- 
pecting a coexistence of a finite condensate fraction and 
broken rotational symmetry (i.e. "supersolidity" ) in 
trapped systems are two- fold as: (i) small trapped sys- 
tems are mesoscopic (ii) inherent incommensurability be- 
tween spherical geometry of the trap and triangular ge- 
ometry of a 2D triangular lattice leads to an effective 
introduction of defects in the lattice, which might lead 
to appearance of superfluidity. To our best knowledge 
quantum systems of two-dimensional trapped dipoles at 
zero temperature so far were not well studied. 

This work is organized as follows. In Section [TT] we 
outline model Hamiltonian and describe methods used in 
simulation, in Section [IIII we describe method of natural 
orbitals used for studying coherence in the system in this 
section we generalized method used in 3D system (see 
[30j . [3l| ) to the case of 2D, in Section HVl we present 
obtained results and in Section [V] we draw conclusions. 

II. MODEL SYSTEM AND NUMERICAL 
APPROACH 

Wc consider a two-dimensional system of N dipolar 
bosons in a harmonic trap with frequency u). Such a 
system is described by the Hamiltonian 

AT 2 N N ,2 

l i i<j L J 

The second term in the Hamiltonian (TTJ) is associated 
with harmonic confinement potential, the third term de- 
scribes the dipolar interaction between particles, m is the 
mass of a particle and ddi P is its dipole moment. It is con- 
venient to use oscillator units in the problem, tha t is to 
measure length in units of oscillator length ciq — yjh/muj 
and energy in units of huj. The dimensionless Hamilto- 
nian is then 

1 N N , 

£=5£(-v?+^)+£^ (2) 

where d = d 2 dip / a^hiu is a dimensionless coupling param- 
eter, which can be interpreted as the ratio of the typical 
energy of the dipolar interaction energy Ei n t = d di /a^ 
and the characteristic energy of a harmonic oscillator 
confinement E tra p = fv-^- 

In order to study system properties we use Monte Carlo 
(MC) methods. A number of MC methods may be em- 
ployed for simulation of quantum systems: variational 
Monte Carlo (VMC), diffusion MC, Path Integral Monte 



Carlo (PIMC), Path Integral Ground State etc. In the 
present work we are interested in ground state (zero tem- 
perature) properties and we do so by means of VMC and 
DMC techniques. 

In VMC method one samples particle distribution ac- 
cording to a chosen trial wave function. By doing that 
it is possible to obtain mean values of an observable by 
averaging its value over the chain of realizations of par- 
ticles coordinates. Method proposed by Metropolis et al. 
[32j is used to generate such a chain. 

Although we do not know the many-body wave func- 
tion exactly, some physically sound ansatz may be used 
for trial wave function so that it depends on particles 
coordinates and some a set of additional parameters. 
Those parameters, referred to as variational parameters, 
are used to minimize energy corresponding to this wave 
function. A variational principle applies, that is the vari- 
ational energy calculated in this way is always larger than 
the ground state energy and equal to it when the trial 
wave function coincides with the ground state wave func- 
tion 

The main idea of the DMC technique is to solve the 
Schrodinger equation in the imaginary time. For suffi- 
ciently long evolution of an arbitrary wave function in the 
imaginary time its projection to the excited states is ex- 
ponentially suppressed compared to its projection to the 
ground state. The problem is that one would rather pre- 
fer to sample the square of the absolute value of the wave 
function (probability density) than the wave function it- 
self, because average values of an observable with diag- 
onal operator is defined as an integral where observable 
is integrated with weight equal to absolute value of the 
wave function squared. But it is impossible to introduce a 
closed real- valued equation defining the time evolution of 
the probability density. A common choice is to calculate 
mixed estimators for observables where observable is in- 
tegrated with weight equal to /(R, t) — ^r-(R)V>o(R-, t), 
with iPt(R) being the trial wave function and ?po(R,t) 
the ground state wave function. A mixed estimator in- 
troduces some bias due to a particular choice of the trial 
wave function, unless the observable commutes with the 
Hamiltonian. An important example of a mixed estima- 
tor being exact is in the evaluation of the ground state 
energy. 

For other observables one can reduce this bias by ex- 
trapolating mixed estimator to the exact one, 

Aextr — ^'A rn j iXe( i Ai r ^ a \. (3) 

Here A m i X 

ed = J iprAipo is mixed estimator of an ob- 
servable A and A tr iai — J V't^/'t is a variational esti- 
mator. This extrapolation is accurate to the second order 
of difference between trial wave function and exact one. 
We applied the extrapolation for all observables that not 
commute with Hamiltonian. 

To reduce the extrapolation error it is important that 
the projection of the trial wave function to the ground 
state is as large as possible. We construct the trial wave 
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function in the Nosanow-Jastrow form 

N N 

* T (R.) = l[f(r i )l[g(r ij ), 

i i<j 



34.95 



(4) 



where R = {r\, r-i, ■ ■ ■ r^} stands for a point in 2N- 
dimensional phase space, f(r) and g(r) are one- and two- 
body correlation terms. The first term in ^t(R-) takes 
into account one-body physics and describes the effects 
of the harmonic oscillator, while the second term intro- 
duces interparticle correlations. For different values of 
the interaction strength d we have used several forms of 
f(r) and g(r). 




A. Gas of dipoles 

An exact expression for the one-body terms /(ft) of 
the wave function @is known in the limit of an ideal 
Bose gas d — > and is given by Gaussians fho(fi) = 
exp (— r^/2a§). We keep this functional dependence form 
for finite values of the interaction strength d: 

f(r) = exp(-ar 2 ), (5) 

with a single variational parameter a, which value is fixed 
by minimizing the variational energy. 

The large distance physics are dominated by the Gaus- 
sian dependence in f(r), so the main requirement for the 
two-body Jastrow term is to describe correctly short- 
range physics. When two particles come close to each 
other the influence of other particles can be neglected and 
the two body Jastrow term in the trial wave function (U) 
can be well approximated by the zero-energy scattering 
solution of the two-body problem, 

g(r) = K (l^dh) ■ (6) 
The resulting trial wave function is, 




FIG. 1. Variational energy as a function of the variational 
parameter fi in a cluster of N = 32 particles with interaction 
strength d = 100. 

used classical Monte Carlo method combined with gradi- 
ent descent optimization to obtain position of sites r| for 
given particle number iV c . After obtaining them we have 
spread classical coordinates with a factor 7 optimized by 
means of variational Monte Carlo. 

Combining Eq. ((5J) with trial wave function of gas 
Eq. ([7]) we obtain the final expression for the crystal trial 
wave function: 




N c N 

xiiE^pH 3 ^- 1 ?) 2 )- ( fl ) 

i i 

The typical dependence of the variational energy on 
parameter /? is shown in Fig[TJ There are two minima. 
The first one at (3 = 0, corresponds to delocalized system 
or a gas state. The second minimum is located at f3 7^ 
and corresponds to a localized system or a crystal. 

C. Gas of strongly correlated dipoles 



B. Crystal of dipoles 

When a crystal is formed the system loses the trans- 
lational symmetry. To take this into account one has to 
use a proper symmetry in the trial wave function. Wc 
introduce a localizing "crystal" term in TWF, 



N c N 

u ( R ) = IlE cx p(^( ri - r J : ) 2 )' 

i i 



(8) 



where (3 is a variational parameter. Each Gaussian term 
describes particle localization close to sites r?, the total 
number of sites in crystal being denoted by N c . We have 



Crystal and gas wave functions provide a good descrip- 
tion in corresponding limits, but their use for intermedi- 
ate values of the interaction strength produce large error 
in mixed estimators. The smaller is difference between 
the trial and ground state and wave functions, the smaller 
is the error. To reduce this difference in region of inter- 
mediate interaction strength, we introduce a shell term 
in the trial wave function, 



s(r) 



n 7 



ai exp (— (r 



M 2 ) 



2 = 1 



(10) 



where r\ is the separation of the i-th shell from the center 
of trap, (Tj is the width of i-th shell and ai is its ampli- 
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tude, 7 is a variational parameter. Parameters of shells 
are optimized so that s(r) (for 7 = 1) is the best approx- 
imation of the ratio of variational and diffusion radial 
distributions. Variational result for the radial distribu- 



tion calculated with s(r) reproduces closely DMC mixed 
estimator calculated without s(r), and it means that the 
trial wave function is better when shell term s(r) is used. 
The resulting trial wave function for this case is 



JV 



v>(R)=n exp( - 



ar,- 



a i exp(-(r-r I s ) 2 / ( T i 2 ) 




(11) 



III. BOSE-EINSTEIN CONDENSATE 

To test the coherence in the system and to calculate 
the condensate fraction we study the spectral properties 
one-body density matrix (OBDM). The OBDM is defined 

as 

p(r,r') = (* t (r')^(r)), (12) 

where ^(r) is the field operator that annihilates a particle 
at the point r. Using single particle states (j>i{r), we 
expand the field operator, 

*(r)=^0 i (r)a i , (13) 

i 

where &i is the bosonic annihilation operator that anni- 
hilates a particle in the state \i). Substituting (fTBf into 
(fT2|) one obtains the spectral decomposition 

p(r,r') = ^*(r) ( (> j (r')N i S ij . (14) 

ij 

The OBDM is diagonal in the single particle state ba- 
sis, and natural orbitals </>i(r) (occupation numbers Ni) 
are its eigenvectors (eigenvalues). The condensate is de- 
scribed by the orbital 4>o(r) with the macroscopic occu- 
pation number and condensate fraction is no = Nq/N. 

The definition of the OBDM operates with eigenfunc- 
tions <P{t) which are functions of four parameters in a 
2D system, so that OBDM cannot be diagonalized using 
standard matrix methods. However the problem for the 
eigenvalues can be simplified if there is a cylindrical sym- 
metry of the problem. In this case the OBDM depends 
on angle 9 = <f>— <f>' between vectors r and r', rather than 
on two angles (j> and (/>' . Furthermore OBDM can be 
expanded in a series of angular momentum components: 

p(r, r') = V Pl (r, r') exp (iW) (15) 

with the projection of the OBDM to the state with an- 
gular momentum I, I = 0, ±1 . . . given by 

pi(r, r') = J d0dr 2 ■ ■ ■ dr N ^*(r,r 2 , . ■ . , r N ) 

xexp(-^0)^(r,r 2) ...,rAr). (16) 



Further, the one-body orbitals can be expanded in terms 
of the angular momentum 

M*) =i?„z(r)$zM, 
H<P) = vfeexp(^), [U) 

where compound index k consists of two indexes l,n. 
Substituting Eq. ^T7\ into Eq. ([T4|) we obtain the repre- 
sentation of the l-th OBDM in terms of natural orbitals 

Mry) = £>«(r)<fe(r')^«- (18) 

i 

Matrices pi(r, r') can be sampled in Monte Carlo simula- 
tion according to Eq. ([TS]) and are "usual" algebraic ma- 
trices understood as functions of two scalar arguments, 
which can be readily diagonalized using standard matrix 
methods. For condensate study we consider only compo- 
nents with angular moment 1 = 0. To solve this equation 
numerically we introduce regularized functions 

u a {r) = 4> a {r)^/r. (19) 

These functions uu(r) are better than <fiu(r) in calcula- 
tions, because they are well behaved near r = 0. In terms 
of u.ii(i), the relation Eq. (|18p reads as 

\frpi(r,r')V? = J2 u U r )Mr')Nu. (20) 

a 

We solve this equation and obtain the condensate wave 
function 

Mr) = uooM/VF. (21) 

IV. RESULTS 

In this section we present results of DMC and VMC 
simulations of the system. In the first part of the section 
we discuss structural properties of the system and in the 
second part we study properties of the condensate. 

A. Structural properties 

We perform a Monte Carlo study of structural and 
energetic properties of dipolar clusters by doing calcula- 
tions with different trial wave functions. We have tested 
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the quality of different trial wave functions for a cluster 
of N — 32 particles in the strongly interacting regime 
d = 200. The results obtained for the energy are sum- 
marized in Table HI 



TABLE I. Ground state energy per particle for N = 32 parti- 
cles for d — 200, energy is measured in units of huj. Number 
in parenthesis show the error on the last digit. 



Method and trial wave function 



Energy per particle 



VMC, gas TWF 

DMC, gas TWF 

VMC, crystal TWF 

DMC, crystal TWF 

VMC, TWF with shell term 

DMC, TWF with shell term 



44.9584(9) 

44.4201(6) 

44.848(1) 

44.4491(5) 

44.8523(7) 

44.4145(8) 




0.000 



FIG. 2. Radial particles distribution and width of the shells 
(inset) for different values of interaction strength d 



One can see that the best variational energy of the cluster 
with strong interactions d = 200 is the one obtained using 
the crystal trial wave function (J5J). Energy of a crystal 
in DMC calculation is larger than both DMC energy of 
a gas and DMC energy of a strongly correlated gas with 
formed shells. The best description of the system is given 
by the wave function shell structured gas (ITT1) . 

DMC calculations of the energy of the cluster with dif- 
ferent interaction strengths showed that the gas phase is 
energetically preferable for values of interaction strength 
smaller than ~ 100, in the region of 100 < d < 400 strong 
correlated gas is energetically preferable and for values 
of d larger than ~ 400 a full crystal is formed. Figure [5] 
shows radial distribution of particles R(r) for different 
values of d obtained with gas trial wave function. One 
finds that R(r) becomes nonmonotonic as the interaction 
increases and a shell structure is formed. We find that 
the shells can be well approximated by a set of Gaussians. 
The width of the Gaussians (see inset in Fig. [2]) satu- 
rates to a constant value when d is large enough. This 
fact can be explained if we assume that dipoles are mov- 
ing around crystal sites in the mean field associated with 
other particles and consider a generalized interaction law 
l/r a , where a > 1. In a homogeneous crystal the first 
not disappearing term of the mean field is quadratic one 



with frequency w 



'mf 



KUr) 



= a(a-l)V(r )/rl 



where is the mean interparticle distance. In the pres- 
ence of the external confinement the size of the cloud is 
proportional to the interparticle distance ro and is such 
that characteristic trapping energy is of the same order 
as the interparticle energy, that is V(tq) ^ uJtrap 

r 2 . Tak- 
ing this into account, one finds that particle displace- 
ment in mean field A is constant for given trap frequency, 
A ~ a(a — l)uJtrap, with a — 3 for dipoles. Results of 
VMC and DMC calculations for radial distribution dif- 
fer for large values of d ~ 100, which means that the 




0.008 



0.006 



0.004 



0.002 



FIG. 3. (x,y) density profile for N = 32 particles in the trap 
for the interaction strength d = 500 



gas trial wave function provides an inadequate descrip- 
tion of the system in the regime of strong interactions. 
We repeated our calculation for trial wave function with 
shell terms and obtained that radial ordering decreases 
the energy, but the lowest energy corresponds to both 
radial and angular localization, i.e. to a crystal phase. 

A contour plot of the density profile of particles is 
shown in Fig. [3] Darker color corresponds to lower den- 
sity of particles. The internal structure of the shells is not 
resolved in Figure [3j as the system can freely rotate due 
to rotational symmetry. In order to extract information 
on the shells ordering we introduce inter- and intra- shell 
order parameter (ctjCHj) with values Ui defined according 
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strength, the condensate is suppressed. The dependence 

T" 




20 40 



60 80 100 

d 



FIG. 4. Diagonal (main graph) and nondiagonal (inset) com- 
ponents of the orientational order parameter (aiOj) with a* 
defined as Eq[22]for cluster with 32 particles 



FIG. 5. Fraction of condensed particles No/N as a function 
of the interaction strength d, number of particles in cluster 
N = 32 



to 

j k 

where <fik is the angle of k — th particle in the j — th shell 
containing Nj particles. For i = j the function (otiCcj) 
is intra-shell order parameter and the larger is its value 
the stronger are internal correlations inside the shell. For 
i =/= j the function (otiOtj) is inter-shell order parameter 
and it measures the correlations strength between shells 
i and j. Results of DMC simulation for a cluster of 32 
particles with a gas trial wave function ([7J are presented 
in Fig. |H One finds that the intra-shell order parameter 
is significantly different from zero for interaction strength 
larger than d ~ 100 and that the inter-shell order param- 
eter is significantly different from zero only for shells 2 
and 3 for interaction strength larger than d ~ 300. Since 
we have used gas trial wave function this correlations are 
not caused by the symmetry of the trial wave function. 

To summarize, we suggest the following scenario of T = 
quantum crystallization. When d is small the gas state 
is energetically favorable, starting from d ~ 10 shells are 
formed in the system, then up to d ~ 300 shells solidify 
and become ordered inside. From d ~ 200 shells start to 
order between each other and around d ~ 400 a complete 
crystal is formed. 



B. Bose-Einstein Condensation and mesoscopic 
supersolid 

We applied the technique of natural orbitals described 
in Section IIII1 to obtain the fraction of Bose-condensed 
particles in a cluster of N = 32 particles with the inter- 
action strength varying from infinitely small {d — 0) to 
intermediately strong (d ~ 100) values. 

To test the method used to estimate the condensate 
fraction we checked that for noninteracting system all 
particles are condensed. At finite values of the interaction 



of the condensate fraction on the interaction strength is 
presented in Fig. [5] Our calculations show that even in 
region of intermediate-strong interactions d ~ 100 the 
fraction of condensed particles significantly differ from 
zero and is order of 40%. At the same time shell structure 
is already formed. Simultaneous coexistence of broken 
rotational symmetry (shells) and a Bose-condensate can 
be referred to as a mesoscopic supersolid. 

V. CONCLUSIONS 

In the present work a two-dimensional system of 
bosonic dipoles in a harmonic trap is studied. Dipoles 
are assumed to be oriented perpendicularly to the two 
dimensional plane and to interact according to repulsive 
dipole-dipole potential. We found that cluster undergoes 
a quantum crystallization for large value of dipole inter- 
action strength (d ~ 400). In the process of zero tem- 
perature transition interaction strength plays the role of 
a control parameter. We found that crystallization con- 
sists of three stages: formation of shells, in-shell ordering 
and ordering between shells. 

We found that quantum Bose-Einstein condensation 
occurs in the system. We investigate dependence of the 
fraction of condensed particles on interaction strength up 
to the intermediate-strong interaction strengths, and find 
that number of particles in condensate can be sufficiently 
different from zero. In this regime spatial ordering (shells 
structure) coexists with Bose condensation, thus a meso- 
scopic supersolid is formed. 
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